Time-periodic solitons in a damped-driven nonlinear Schrodinger equation 
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Time-periodic solitons of the parametrically driven damped nonlinear Schrodinger equation are 
obtained as solutions of the boundary-value problem on a two-dimensional spatiotemporal domain. 
We follow the transformation of the periodic solitons as the strength of the driver is varied. The 
resulting bifurcation diagrams provide a natural explanation for the overall form and details of the 
attractor chart compiled previously via direct numerical simulations. In particular, the diagrams 
confirm the occurrence of the period-doubling transition to temporal chaos for small values of dis- 
sipation and the absence of such transitions for larger dampings. This difference in the soliton's 
response to the increasing driving strength can be traced to the difference in the radiation frequencies 
in the two cases. Finally, we relate the soliton's temporal chaos to the homoclinic bifurcation. 
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I. INTRODUCTION 

The application of a resonant driving force is an effi- 
cient way of compensating dissipative losses in a soliton- 
bearing system. If the dissipation coefficient and driving 
strength are weak, and the driving frequency is just below 
the phonon band, the amplitude of the arising oscillating 
soliton is governed by the nonlinear Schrodinger equation 
with damping and driving terms. 

The damped-driven nonlinear Schrodinger equations 
exhibit localised solutions with a variety of temporal be- 
haviours, from stationary to periodic and chaotic. There 
is a whole range of analytical and numerical approaches 
to the study of stationary and steadily travelling solitary 
waves. As for the solitons with nontrivial time depen- 
dence, such as periodic, the direct numerical simulation 
has remained an exclusive means of obtaining these so- 
lutions and classifying their stability. 

The shortcoming of this method is that simulations 
capture only stable regimes. This means that the actual 
mechanisms and details of the transformations of solitons 
(which are bifurcations involving both stable and unsta- 
ble solutions) remain inaccessible. Neither can simula- 
tions be used to identify coexisting attractors in cases of 
bi- or multistability. 
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In this paper we pursue a different approach to the 
analysis of these hidden mechanisms. Instead of direct 
numerical simulations, the time-periodic solitons are de- 
termined as solutions of a boundary-value problem for- 
mulated on a two-dimensional spatio-temporal domain. 
The advantage of this approach is that it is potentially 
capable of furnishing all solutions — all stable and all 
unstable. 

The particular equation that we are concerned with 
here, is the parametrically driven damped nonlinear 
Schrodinger, 

iipt + i^xx +'2\ip\'^tp - ^ = hip* - i-f>jj. (1) 

In Eq. ([!]), 7 > is the damping coefficient and h the 
amplitude of the parametric driver, which can also be 
assumed positive. Equation ([T]) was used to model a vari- 
ety of resonant phenomena in nonlinear dispersive media, 
including the nonlinear Faraday resonance in a vertically 
oscillating water trough [T}^3|, formation of oscillons in 
granular materials and suspensions [4], synchronisation 
in parametrically excited pendula arrays |51 [B], phase- 
sensitive amplification of light pulses in optical fibers [7] 
and propagation of magnetization waves in an easy-plane 
ferromagnet placed in a microwave field IHl [9]. The 
same equation governs the amplitude of breathers in a 
variety of systems reducible to the parametrically driven 
damped sine-Gordon |8J and the (j>'^ [T^ equation. (For 
more contexts, see [TT].) 

The studies of localised states in these mechanical, hy- 
drodynamical, magnetic and optical systems have been 
focussing on structures oscillating periodically, i.e., with- 
out incommensurable frequencies in their spectrum. The 
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quasiperiodically oscillating states (e.g. states whose fun- 
damental oscillation is modulated by smaller frequencies) 
have been set aside as complex and atypical. The pe- 
riodically oscillating states in the physical systems re- 
ducible to equation ([T]) are described by the stationary 
solutions of this equation. Accordingly, all previous anal- 
yses of localised solutions of Eq.Q have been confined 
to its stationary solitons. One of the conclusions of the 
present project, however, will be that the time-periodic 
solitons of Eq.Q populate a significant part of its at- 
tractor chart. Therefore, the corresponding quasiperi- 
odic localised states should play a much bigger role in all 
physical applications of Eq.Q than it was assumed so 
far. 

The aim of the present work is to follow the transfor- 
mations of temporally periodic solitons of equation ([T]) 
as its parameters are varied, identify the arising bifur- 
cations and eventually verify and explain the attractor 
chart for this equation which was compiled using direct 
numerical simulations in Ref. [H]. In the second part of 
this project [11] we will complement this one-soliton chart 
with a chart of two-soliton attractors. 

The paper is organised as follows. Section |ll] (which 
continues this introduction) contains background infor- 
mation on stationary solitons and their transformations 
as the parameters of the equation are varied. It is these 
transformations that we will be verifying and studying 
further in the subsequent sections. 

Next, sections |III| and [TV] present mathematical tech- 
niques we employ in this project: in section [lIT| we outline 
our method of obtaining periodic solitons whereas section 
IV describes our approach to the analysis of their stabil- 
ity. Section |VI| introduces a theoretical framework for 
the treatment of radiation from the oscillating soliton. 

The central results of this study are obtained using nu- 
merical methods; these are reported in section|V] Here we 
present bifurcation diagrams for the time-periodic free- 
standing soliton in various damping regimes. 



where 



The paper is concluded by section VII where the results 



of the numerical study are discussed and interpreted. 



II. STATIONARY AND OSCILLATORY 
SOLITONS: THE BACKGROUND 

Localised stationary or periodic solutions of Eq. ([T]) ex- 
ist only if /i > 7. When h > /icont (7), where 

^cont = + 7^1 

any localised solution is unstable to continuous-spectrum 
perturbations. The evolution of this instability leads to 
spatiotemporal chaos. 

Two stationary soliton solutions of Eq.Q are well 
known. One soliton (denoted in what follows) ex- 
ists in the parameter range j < h < /icont (7) and has the 
form 
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This solution is unstable for all h and 7 [2 |8j. We 
are mentioning this unstable object here because it will 
reappear below as a constituent in stationary multisoli- 
ton bound states. We will also be recalling this soliton 
when interpreting complex temporal behaviour of time- 
periodic solutions. 

The other stationary soliton exists for all h > we 
denote it ip^: 



ip+{x) — A+ exp(— i0-|_) sech(A_|_a;), 



where 
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tp-{x) — j4_ exp(— i0_) aech{A^x), 



(2a) 



The stability properties of this soliton depend on 7 and 
h [S. When 7 > 0.356, the ip+ soliton is stable for all 
h in the range j < h < /icont(7)- When 7 < 0.356, on 
the other hand, the soliton ([3| is only stable for 7 < ft, < 
/iHopf(7), where the value /iHopf(7) lies between 7 and 
/icont(7) (see the curve labelled 1 in Figjl]). As we increase 
h past ftHopt(7) keeping 7 fixed, the stationary soliton 
loses its stability to a time-periodic soliton [51 [T3] . The 
transformation scenario arising as h is increased further 
depends on the choice of the (fixed) value of 7. 

The numerical simulations ''12] (also [14J) indicate that 
for 7 smaller than approximately 0.25, the periodic soli- 
ton follows a period-doubling route to temporal chaos. In 
a wide region of h values above the chaotic domain, the 
equation does not support any stable spatially-localised 
solutions. In this "desert" region, the only attractor 
found in direct numerical simulations was the trivial one, 
-0 = 0. Finally, for even larger values of h, the unsta- 
ble soliton seeds the spatio-temporal chaos [12] (see also 

m)- 

As h is increased for the fixed 7 greater than approx- 
imately 0.275, the soliton follows a different transforma- 
tion scenario. Here, the period-doubling cascade does 
not arise and the soliton death does not occur. The pe- 
riodic soliton remains stable until it yields directly to a 
spatio-temporal chaotic state [T2] . 

In a short intermediate range of 7-values, 0.25 < 7 < 
0.275, we have a combination of the above two scenarios. 
The increase of h for the fixed 7 results in the period- 
doubling of the soliton, culminating in the temporal- 
chaotic regime, which is followed by the soliton death. 
As we continue to raise h, an inverse sequence of bifurca- 
tions is observed which brings the stable single-periodic 
soliton back. On further increase of h, it loses its stability 
to a spatio-temporal chaotic state [T^ . 

The attractors arising for various h and 7 values are il- 
lustrated by Figjl] which we reproduce from Ref . [12J . We 
should emphasise that this attractor chart has been com- 
piled using direct numerical simulations of equation ([I]), 
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spatio-temporal chaos 




FIG. 1. The single-soliton attractor chart for equation ([I]) 
compiled by direct numerical simulations 12 . Stationary 
solitons serve as attractors in the blank area above the h — ^ 
line. The area where only the trivial attractor is available is 
marked by the empty triangles. Empty circles indicate sta- 
ble periodic solitons; the black triangles label spatio-temporal 
chaos. Other symbols mark stable higher-periodic and tem- 
porally chaotic solitons. For details see [T^ . 

with a particular choice of initial conditions. (The initial 
condition was chosen in the form of the unstable soliton 
— in the parameter range where it is unstable.) It 
is an open question, therefore, how robust this chart is. 
Would it change if simulations started with a different 
initial condition, or if one used a numerical scheme with 
different parameters (such as the spatial interval length, 
the number of the Fourier modes, the full time of simula- 
tion, temporal stepsize, etc)? In particular, would qual- 
itative features of the chart survive, e.g. the coexistence 
of two transformation scenarios, appearance of the desert 
region, and the peculiar shape of the periodic-attractor 
domain? One of the aims of the present work is to answer 
these questions. 

III. PERIODIC SOLITONS AS SOLUTIONS OF 
A BOUNDARY- VALUE PROBLEM IN 2D 

Instead of solving equation ([l]) with some initial condi- 
tion and determining the resulting attractor by running 



the computation for a sufficiently long time, we were 
searching for periodic solutions by solving Eq.Q as a 
boundary-value problem on a two-dimensional domain 
(—00,00) X (0,T). The boundary conditions were set as 

ip{x,t)^0 as x ±00, (4) 

and 

^{x,t + T) ^^{x,t). (5) 

The value of T was not available beforehand; the period 
was regarded as an unknown, together with the solution 
ip{x,t). 

The periodic solutions were continued (path-followed) 
in h for the fixed 7. We employed a predictor-corrector 
continuation algorithm [T5] with a fourth-order Newto- 
nian iteration at each h. A finite-difference discretization 
with the stepsize Aa: = 0.05 was used on the interval 
i-L,L) = (-50,50). 



IV. STABILITY OF PERIODIC SOLUTIONS 

A. Floquet multipliers 

Let ipoix, t) = TZ{x, t) + iX{x, t) be a spatially locaHsed, 
time-periodic solution. Letting ip{x,t) = ■0o(x,t) -I- 
u{x,t) + iv{x,t) and linearising ([T]) in the small pertur- 
bation u + iv, we obtain 

Jwt = {H- jJ)^, (6) 

where w = w(a;,t) is a two-component column-vector 

J is a skew-symmetric matrix 

and "H is a hermitian matrix-differential operator 



n 
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-mi 



-mi 
i-h-2n^ - 61^ 



(7) 



The solution to Eq.([6]) with an initial condition w(a;, 0) 
can be written, formally, as w{x,t) = A4t'^(x,0), where 
the evolution operator Ait acts on (vector-) functions of 
X but depends, parametrically, on t. A fundamental role 



I 

is played by eigenvalues of the monodromy operator Mt- 

Mryix) = fiy{x). (8) 

Here is the evolution operator Mt evaluated at 

t = T. The eigenvalues fi are usually referred to as the 
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Floquet multipliers and the exponents A, where fi 
as the Floquet exponents. According to the Floquet the- 
ory, for each A there is a solution w{x,t) such that 



w{x,t) = e^*^p{x,t), 



(9) 



where p{x,t) is periodic with the period of the solution 
tpQ-. p{x, t + T) = p{x, t) for all t. 

It is useful to establish the relation between the 
evolution operator Mt and symplectic maps. Letting 
■w{x,t) = e~''*\v{x,t), Eq.Q is cast in the form 



(10) 



Since T-l is hermitian, equation ( 10 ) is a hamiltonian sys- 



tem (with a quadratic Hamilton functional). The solu- 
tion to ( 10 1 with an initial condition w(a;, 0) is, therefore. 



w(a:, t) — 5tw(x, 0), 

where St is a symplectic map. Thus the evolution op- 
erator Mt and, in particular, the monodromy operator 
Mt, are related to symplectic maps: Mt = e~^*St and 



Ml 



(11) 



We will use this relation below, in order to explain sym- 
metries of the set of eigenvalues of the operator Mt- 



B. Spectrum structure 

Stability properties of stationary solutions are deter- 
mined by eigenvalues A of the operator {J~^'H — 7). 
The corresponding Floquet multipliers are given by e^^ 
(where T can be chosen arbitrarily for a stationary solu- 
tion). Therefore the stability eigenvalues A are nothing 
but the Floquet exponents. Each of the two stationary 
solitons '0+ and ip- has just one zero stability eigenvalue 
(i.e. just one exponent A = 0). This zero mode origi- 
nates from the translation invariance of equation ([l]) — 
this is the only continuous symmetry the damped-driven 
equation has. 

At the point of the Hopf bifurcation of a station- 
ary solution (a free-standing soliton -0+ or a stationary 
bound state of solitons), two complex-conjugate eigenval- 
ues A, A* cross into the ReA > half-plane. Irrespectively 
of what one takes for the value of T in this case, the cor- 
responding Floquet multipliers — e^^ and ^* — 
cross through the unit circle on the complex //-plane. If 
we, however, take T to be equal to the period of the 
periodic solution bifurcating off at this point, i.e., let 
T ~ 27r/ImA, the Floquet multipliers cross through the 
unit circle exactly at /i = 1. Adding the unit multiplier 
associated with the translation invariance, we conclude 
that any stationary solution has 3 unit Floquet multipli- 
ers at the point of its Hopf bifurcation and, accordingly, 
the detaching periodic solutions should also have 3 unit 
multipliers at that point. 



As we continue the periodic solution ipoix-, t) away from 
the point where it was born, two of the unit Floquet mul- 
tipliers persist (while the third one moves away along the 
real axis) . These two unit multipliers are associated with 
the translation invariance in x and periodicity in time, 
respectively. Indeed, if we substitute ^pQ{x,t) back in 
Eq.([T]) and differentiate the resulting identity with re- 
spect to X and t, we will obtain Eq.Q with w = dxipo 
and w = dfipo, respectively, where the two-component 
vector 



n{x,t) 
I{x,t) 



This means that equation Q has two periodic solutions 
and the monodromy operator Q has two eigenvalues 
/i = 1 with eigenfunctions w = o^V'o and w = dtipa, re- 
spectively. These two unit eigenvalues will routinely arise 
in our stability analysis of periodic solutions of Eq.Q. 

It is not difficult to show that real eigenvalues of Mt 
will always arise in pairs whereas the complex eigenvalues 
will always come in quadruplets. Indeed, the monodromy 
operator is related to a symplectic map by the relation 
(111. Real eigenvalues of symplectic maps are known to 



come in pairs; if is a real eigenvalue of St, then so is 
v^^. Complex eigenvalues of St come in quadruplets; if 



f is such an eigenvalue, then so are i' 



and (i/* 



|16) . The relation (11) implies then that if is a real 



eigenvalue of the evolution operator Mt, then so is p,, 
its inverse with respect to the circle of radius e~'^'^: 



-27T 



In a similar way, if /i is a complex eigenvalue of A^^, then 
/i*, fl, and fl* are eigenvalues as well. 

We will be referring to fl simply as the mirror image 
of fl. Out of the two Floquet multipliers fi and /t only 
one can cross the unit circle whereas its mirror image 
will always be confined inside a smaller circle of the ra- 
dius e"'''^ . Consequently, when analysing the motion of 
multipliers on the complex plane and resulting stability 
changes, we will be focussing on /i with moduli greater 
than e~'^'^ and ignoring those with < e""*""^. 

Finally, we discuss the location of the continuous spec- 
trum. As |a;| — > 00, the operator (J~^'H — 7) becomes 
a matrix differential operator with constant coefhcients 
whose spectrum is easily determined. Namely, when < 
1, the spectrum consists of all A of the form A = —7 ± ip, 
with y/l — h'^ < p < 00. The corresponding Floquet 
multipliers fill in a circle of the radius e~'''^ on the com- 
plex /x-plane. When > 1, the continuous spectrum of 
{J~^H — 7) fills the vertical line Re A = —7 and an inter- 
val on the real axis, —7 — V/i^ — 1 <A< —7 + ^/h"^ — 1. 
The corresponding Floquet multipliers fill the circle of 
the radius e~'''^ and in addition, an interval on the real 
axis: e-(^+^'^)^ < < e-(7-v^)T_ 
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FIG. 2. The period of the periodic solutions with 7 = 0.30 
and 7 = 0.35. The sohd curves show stable and the dashed 
ones unstable branches. 



C. Numerical stability analysis: the method 

To find out whether equation (|6| admits solutions of 
the form (|9| with ReA > 0, we expand u(x^ t) and v{x, t) 
in the Fourier series on the interval (— i, L). The bulk of 
our eigenvalue calculations were done with N = 100 but 
we went up to iV = 250 when the eigenf unctions have 
shown variations on a small scale. 



V. NUMERICAL STUDY 
A. Strong damping: 7 = 0.30 and 7 = 0.35 

We explored 7 = 0.30 and 7 = 0.35 as two rep- 
resentative sections of the attractor chart in its right- 
hand part, where numerical simulations had detected no 
period-doubling bifurcations. The transformation of the 
solution as it is continued in h is similar in the two cases; 
see Fig(2| 

The top-end point of each of the two curves in figure [2] 
corresponds to the stationary single-soliton solution ipj^.. 
The underlying value of h equals 0.385 for 7 = 0.30 and 
0.7500 for 7 = 0.35. At these /i, the stationary tpj^ soli- 
ton undergoes a Hopf bifurcation and a stable periodic 
solution is born. At the starting point of each curve, 
the spectrum of the periodic soliton includes three unit 
multipliers /xi_2.3 = 1- As /i is increased, the eigenvalue 
/i3 moves inside the unit circle whereas /ii_2 remain at 1. 
Meanwhile, five real positive eigenvalues iJ,i^...fi detach, 
one after another, from the continuous spectrum. One 
of these (/ig) later returns to the continuum while the 
other four eigenvalues move towards the unit circle. At 
some point the eigenvalue ^4 collides with /i3 producing 
a complex pair which, however, later returns to the pos- 




FIG. 3. (Color online) The absolute value of the periodic 
solution with large 7. In this plot, 7 = 0.35. "Motion pic- 
tures" (a) and (c) have been taken at the same value of h but 
correspond to the different branches of the diagram Fig[2] (a) 
h = 0..85, r = 2.45 (top branch) and (c) h = 0.85, T = 2.34 
(bottom branch). The evolution (b) corresponds to the turn- 
ing point: h — 1.0185, T = 2.21. In each case several periods 
of oscillation are shown. 

itive real axis. At the turning point hsn, two positive 
eigenvalues, fis and /i4, cross through the unit circle (al- 
most simultaneously). This is where the periodic solution 
loses its stability. Numerically, the turning-point value is 
hsn = 0.8761 for 7 = 0.30 and h^n = 1.0186 for 7 = 0.35. 
On the unstable branch, the spectrum includes two pos- 
itive eigenvalues /Lt3^4 > 1, two unit eigenvalues ^11,2 — 1, 
and one positive eigenvalue close to unity {p^ < 1). In 
addition, two more positive eigenvalues approach the unit 
circle from inside as we continue away from the turning 
point. 



It is fitting to note here that the two turning point 
values, hsn — 0.8761 and /ign — 1.0186, are in a good 
agreement with the boundaries of the periodic-attractor 
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existence domain established previously. Namely, the di- 
rect numerical simulation of Eq.([T]) gave values close to 
0.86 and 1.01 for 7 = 0.30 and 0.35, respectively [H]. 

The end point of the dashed curve {h = 0.61 for 
7 = 0.30 and h — 0.760 for 7 — 0.35) corresponds to a 
stationary complex of solitons. In the 7 = 0.35 case, this 
solution has three separate humps in its real and imagi- 
nary part. Comparing (the real and imaginary parts of) 
each of the three humps to (the real and imaginary parts 
of) the free-standing solitons V'+ and V'-; we conclude 
that the complex consists of the ip^ soliton in the middle 
and two solitons V'- on its sides; that is, the stationary 

complex should be interpreted as ^( i ). The spectrum 

of discrete eigenvalues of this stationary complex is close 
to a union of the eigenvalues of the soliton ipj^. and those 
of two solitons ijj- . Namely, we have three unit eigenval- 
ues fJ.1.2,5 = 1 (resulting from the translation modes of 
two ip^^s and one ip^) and two eigenvalues /^6,7 close to 
unity (contributed by the ip+ soliton which is close to its 
Hopf bifurcation point). In addition, the two "0- solitons 
contribute two real positive eigenvalues //3_4 > 1. 

In the case of 7 = 0.30, the two side peaks in the real 
part of the end-point solution are seen to have merged 
with the central peak (whose height coincides with the 
height of the soliton in its real part). On the other 
hand, the central peak in the imaginary part of the so- 
lution is seen to have disappeared while the two lateral 
peaks have the same height as the soliton ip- in its imag- 
inary part. This indicates that the solution can still be 

interpreted as the ■;/'( i ) complex, albeit a tightly bound 

one. 

Representative solutions are shown in figure [3j Near 
the leftmost point of the curve pertaining to 7 = 0.35 in 
Figj2] the periodic solution looks like a single soliton with 
a periodically oscillating amplitude and width [Fig Sjja)]. 
The soliton is emitting radiation waves; however, the ra- 
diation is decaying rapidly as |x| — ?> cxi. As we move 
further along this curve in Figj2] the amplitude of os- 
cillations as well as the intensity of radiation increases 
[Figj3]jb)]. The shape of the oscillating solution evolves 
into a three-hump structure. Near the end point of the 
curve, the amplitude of oscillations decreases [Figj3](c)] 
and we arrive at the stationary three-soliton complex. 

Thus, periodic solitons with large 7 connect (in the 
sense of paths in the parameter space) stationary solitons 
and their complexes, with the connection points provided 
by the Hopf bifurcations of the latter. It is appropriate 
to mention a recent publication |17j where a similar or- 
ganisation of the solution manifold was reported for the 
Benjamin-Ono equation. 



B. Weak damping: 7 — 0.1 and 7 = 0.2 

In the weak damping regime we explored two repre- 
sentative values of 7, 7 = 0.10 and 7 = 0.20. The two 
cases exhibit similar bifurcation diagrams (Fig. [4]) which 
are, however, very different from the strong-damping di- 
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FIG. 4. The period of the periodic solutions with 7 = 0.10 
and 7 = 0.20. The solid curves show stable and the dashed 
one unstable branches. 



agrams in Fig [2] 

As before, the starting point (the left-end point) of 
each branch corresponds to the stationary ■0+ soliton. 
At h = 0.1250 (for 7 = 0.10) and h = 0.2275 (for 
7 = 0.20) the stationary soliton undergoes a Hopf bi- 
furcation and a periodic solution is born. In Fig[4| this 
solution is marked Tl. The spectrum of the Tl solution 
at the starting point of the curve has three unit mul- 
tipliers, A*i,2,3 = l- As/iis increased, fii moves along 
the real axis inside the unit circle and immerses in the 
continuous spectrum at /i = e~''''^ . For even a larger h, 
a pair of complex-conjugate multipliers with real parts 
close to —e~'^'^ detaches from the continuum before con- 
verging on the negative real axis. One of the resulting 
real eigenvalues moves towards — e~^^ and rejoins the 
continuum. The other one moves away from the origin 
and eventually passes through —1. This is a signature of 
the period-doubling bifurcation: the original solution of 
period T loses its stability and a stable solution with a 
period 2T is born. (This happens at h2T = 0.1463 and 
= 0.2663 for 7 = 0.10 and 7 = 0.20, respectively.) 

The spectrum of the newly born T2 solution evolves 
similarly to the spectrum of Tl. It is also born with three 
discrete multipliers /ii,2,3 — 1- As /i increases, fj,i moves 
along the real axis; immerses in the continuous spectrum; 
reappears at —e~'^'^, and eventually passes through —1 at 
some h = h4T- (The value /i4t equals 0.1511 for 7 = 0.10 
and 0.2725 for 7 = 0.20.) At this point a new stable 
branch is born, with the period equal to two periods of 
T2. This new, period-four branch is not shown in Fig. [4] 

As for the (unstable) Tl solution, the corresponding 
T{h) curve is snaking up, with the period growing with- 
out bound. The solution approaches a homoclinic orbit, 
connecting the stationary -!/;_ soliton to itself. In the next 
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subsection, we will discuss this phenomenon in more de- 
tail. 

We note that the bifurcation diagrams of Fig. |4] are 
in agreement with direct numerical simulations of Eq. ([l]) 
with 7 = 0.10 and 7 = 0.20 [12j that revealed the period- 
doubling transition to temporal chaos in soliton's dynam- 
ics. It is therefore natural to expect that our TA branch 
will also undergo a period-doubling bifurcation after a 
short interval of stability, and similarly the branches that 
it births would undergo a sequence of period-doubling bi- 
furcations. We did not have computational capacity to 
verify this using our two-dimensional continuation ap- 
proach. It is worth mentioning, however, that the values 
of h at which higher-periodic and temporally-chaotic at- 
tractors were observed in simulations (0.16 and 0.28 for 
7 = 0.1 and 0.2, respectively [T^]) are close to the /i4T 
given above. 

Finally, we have not been able to perform an accurate 
numerical continuation of solutions with 7 = 0.05 and 
smaller. In this case the solution consists of a finite- 
extent soliton riding on a (second-harmonic) oscillatory 
background of a non-negligible amplitude, which shows 
only a very slow spatial decay as |x| — ?> 00. [Figjsjjb) gives 
an idea of the shape of the solution in that case, for most 
h.] In order to obtain this solution under the boundary 
conditions ip{ztL,t) = 0, one has to enlarge the length of 
the spatial side of the domain of computation, {—L, L) x 
(0, T). This quickly saturates the computational capacity 
available. 



C. Intermediate damping: 7 — 0.265. 

Finally, we investigate the transformation of the pe- 
riodic solutions with 7 lying between the strong- and 
weak-damping ranges. As a representative value of such 
"borderline" damping, we take 7 = 0.265. According 
to the numerical simulations, the direct period-doubling 
cascade is followed by an inverse cascade here, resulting 
in a peculiar shape of the periodicity region on the (7, h) 
plane (Figjl]). Our aim is to provide an explanation for 
this phenomenon on the basis of the transformations of 
the periodic solutions as h is continuously varied. 

The results of our numerical continuation are sum- 
marised in Figj5] The bifurcation diagram consists of 
three branches, the second and third of which arise as a 
result of the period-doubling bifurcation of the solution 
on the first branch. The transformation of the solution 
as we move along this period-one branch [the bottom 
branch in Figjs] denoted Tl] is similar to the transfor- 
mation of the solution along the curves shown in Figj2] 
The starting point of the curve (its leftmost point) corre- 
sponds to the stationary solution ip^; the end point cor- 
responds to the stationary complex ■!/)(- ^ ). The major 

difference from the case of large 7 (7 = 0.3 and 7 = 0.35) 
occurs when a real eigenvalue crosses through — 1 as /i 
is increased through /i = 0.44. As /i is increased past 
0.44, this negative eigenvalue continues to grow in abso- 
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FIG. 5. The period of the periodic solution for 7 = 
0.265. Solid curves show stable and the dashed ones unstable 
branches. 



lute value, reaches a maximum, and then starts to de- 
crease. As h is increased through h = 0.75, the negative 
eigenvalue crosses through —1 once again, this time in 
the direction of decreasing modulus. The described be- 
haviour of the negative eigenvalue corresponds to period- 
doubling bifurcations at h = 0.44 and h = 0.75, where 
periodic solutions with double period are detached. For 
larger h the solution is stable — all the way to the turn- 
ing point, where two real eigenvalues cross out through 
the unit circle. 

The spectrum of the first double-periodic solution (de- 
taching at h = 0.44 and denoted T2 in Figj5]) includes 
two unit eigenvalues, /ii_2 = 1- When h — 0.44, there is 
also a negative eigenvalue = — 1. As /i is increased 
from 0.44, ^3 moves inside the unit circle, then reverses 
and crosses through —1 once again (at h^T = 0.495). At 
this point, a new periodic solution T4 is born, with the 
period equal to double the period of the solution T2 — 
roughly four times the period of Tl. (This "period- four" 
solution is not shown in Figjsj) The value h^T = 0.495 
coincides with the largest value of h at which higher- 
periodic attractors were seen in simulations |12j . 

Returning to the T2 solution, it has two positive eigen- 
values //4^5 in addition to the eigenvalues /ii,2,3- As the 
T2 curve in Figj5] is traced around the turning point at 
h = 0.543, these two eigenvalues cross, almost simulta- 
neously, through the unit circle so that /^4^5 > 1 on the 
upper branch of the curve. As we continue further along 
the upper branch, the period of the solution grows: the 
solution develops a long epoch where it remains very close 
to the stationary soliton Each period now consists 
of two phases: the solution performs a rapid oscillation 
with its spatio-temporal profile close to that of the Tl so- 
lution, followed by a slow passage through the bottleneck 
near the ip- soliton [see Figj6]. 



8 




1.5 Re\|/ 



FIG. 6. (Color online) (a) The absolute value of the double- 
periodic solution found on the left upper branch in Fig[5] 
Here 7 = 0.265, h = 0.502 and T = 22.985. A rapid os- 
cillation is followed by a long quasistationary epoch where 
the solution remains close to the tp- soliton. (b) The phase 
portrait of this solution taken at a; = 0. Shown is lmtp{0,t) 
vs Ketf}{0,t). The filled and open circle represent two fixed 
points: the stationary solitons tp+ and tp-- 



The second period-2 solution (denoted T2 in Fig|5|, 
which detaches at h — 0.75, remains stable for 0.731 < 
h < 0.75. As with the first double-periodic solution T2, 
the period of T2 grows as we move along the branch. The 
solution changes similarly to what we have described in 
the previous paragraph: a rapid oscillation is followed by 
a long quasistationary epoch when the solution is very 
close to the stationary soliton ip^. 

Thus, periodic solutions on the T2 and T2 branches 
approach a homoclinic orbit: an infinite-period solution 
which tends to the stationary soliton ■(/;_ as t — )■ —00 
and i — > cx). The same transformation scenario was de- 
tected in the weak-damping case (7 < 0.25) where the 
Tl (and, presumably, the T2) branch was seen to snake 
up to T = 00 (FigQ. At the point h — h^o where the 
period becomes infinite, the periodic solution undergoes 
a homocilinic bifurcation which may serve as the source 



VI. RADIATION FROM THE OSCILLATING 
SOLITON 

A. The long-range radiation indicator 

Like their stationary counterparts, the oscillatory soli- 
tons strike the balance between the energy fed by the 
driver and the energy lost to the dissipation. This can 
be quantified, for example, by considering the integral 
N — IV^pda;. [When the nonlinear Schrodinger equa- 
tion is employed to describe planar stationary waveg- 
uides, this integral measures the total power of light cap- 
tured by the {—R, R) section of the waveguide. There 
are energy- related interpretations of this integral in other 
contexts as well.] 

Equation ([T]) yields 



N ^2h I IV'I' 



sm{2e)dx - 2-fN + $1 



(12) 



where we have decomposed tp as \ilj\e *^ and defined 



(13) 



The first term in the right-hand side of ( 12 ) gives the rate 



at which the energy is pumped into the soliton while the 
second one quantifies the damping rate. Assuming that 
the interval {—R, R) is large enough to contain the core of 
the soliton, the last term in ( 13 ) measures the radiation 
fiux through its endpoints. 

The radiation through the points x = ±i? outside the 
soliton's core is governed by the linearisation of ([I]) whose 
dispersion relation is 



{lJ - i^f = (fc2 + 1)2 _ e 



(14) 



The radiation waves emanate from the oscillating soliton 
which plays the role of a pacemaker for these waves. The 
pacemaker has a a period T and will sustain waves with 
frequencies a; = 2im/T {n — 1,2,...); hence uj should be 
taken real in ( 14 ). For real w, equation ( 14 ) has two pairs 
of complex roots k — ±(j3i.2 + iqi,2), with the imaginary 
components 



91,2 



/l S 1 r- 
\ -±- + 
V 2 2 2^ 



2(Q±5), 



(15) 



where 



and 2V - 7^ -l-tj^. 

The decay rate qi is always greater than q2 and 
bounded from below {qi > 1); it accounts for rapidly 
decaying core of the soliton. Outside the core, the decay 
rate crosses over to 92, and if (72 is small the solution en- 
ters the long oscillatory wing decaying in proportion to 



-92^1 



Periodic solutions that we have reported in this 



of chaos (see section VII A below) 



paper all have ui ^ 1; for these, 
small only if 7 is small. 



the exponent 52 can be 
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Assuming that 7 is small, the expression ( 15 1 simpli- 
fies. Defining the quantity 



(16) 



we have, in particular, 



92 



as^^O (.>0); (17a) 
(17b) 



92 



^1/2 as - ^ 0; 

1 + cr 7 



92 ^ (-cr)^/-^ as - ^ (cr < 0). (17c) 
cr 

Thus 92 is small if 7 is small while cr is positive. In 
particular, if cr is much greater than 7, the exponent 
92 ^ l/V'^ while if cr is much smaller than 7, we have 
92 ~ y/j- The decay rate stays small (^ -^Z— cr) when cr 
becomes negative — as long as it remains much smaller 
than 7 in absolute value. However, as a grows to larger 
negative values, the exponent 52 grows as well. Hence 
the quantity cr(tL') serves as an indicator of whether the 
long-range radiation with the frequency w can be excited 
(i.e. whether (72 is small) or not. 



B. Radiation frequency selection 

The harmonic-wave solution of the linearised equation 
with the decay rate (72 is 



4> — rj [ cos (f) + i 



a sin (j) — l3 cos 



(18) 



where the phase <f> — ut — p2X, the wavenumber p2 — 
— T/(2(72), and the coefficients 

a = uj{S + K)+ 7r, /3 = 7(5 + /i) - cjT, 

and T — \JQ,~V ■ The amplitude r\ = rjioj) is arbitrary 
as far as the linearised equation is concerned but acquires 
a specific value when ( 18 ) is used to represent radiation 
from the soliton. 

If R is large enough, the solution ^j{x,t) with x near 
±R is described by ([I8]). Substituting ([l8| in ([l3]) we 
obtain the flux through the points x — ±R: 



^±R)=T\v\ 



, T uj{S + h)+ -fT _ 

q2 {S + hY+T^ ^ 



292 -R 



(19) 



Assume, first, that (72 is of order 1; this occurs, in par- 
ticular, when 7^1. Since R is chosen to be large, the 
fiux through the points x — ±R in this case is exponen- 
tially small: strongly damped radiation waves decay so 
quickly that they do not reach beyond the core of the 
soliton (and essentially form part of it). The exponen- 
tial factor in ( 19 ) is not exponentially small only if q2 is 



small; this, in turn, may only happen when the damping 
is weak. 



These considerations have a simple physical interpre- 
tation. When 7 is large, the soliton dissipates most of 
its excess energy within its core; this is accounted for by 
the second term in (12). On the other hand, when 7 



is small, the dissipation within the core is insufficient to 
balance the first term in (12 1. The long-range radiation 



may provide an alternative energy-draining mechanism 
in this case, but this is only available when (72 is small 
which, in turn, is indicated by cr > 0. 

If the period of the solution is small enough for the 
quantity cr(27r/T) to be positive and much greater than 
7 (which is assumed small), the decay rate q2 is close 
to zero. Therefore, the flux ( 19 1 is not exponentially 



small and the first-harmonic radiation is not blocked in 
this case. As T is increased so that — j'^ + (27r/T)^ 
drops below 1, the decay rate (72 grows to values of order 
1. According to eq.(19), the first-harmonic radiation will 
be mostly suppressed in this case (i.e. damped outside 
a small neighbourhood of the soliton's core). However, 
a{uj) with uj = An/T will be still positive; this means 
that the tails will be dominated by the second-harmonic 
radiation now. Once the second harmonic is suppressed, 
the third harmonic will take over, and so on. 

We note that the second-harmonic radiation is present 
not only when the first-harmonic radiation is suppressed. 
However, the amplitude rj{2nn/T) typically scales as e" 
and so Eq.(19) implies that the second-harmonic fiux is 



much weaker than the first-harmonic fiux in situations 
where both channels are open. 



C. First-harmonic radiation as a stabilising agent 

We routinely calculated the indicator cr(27r/T) as we 
continued our periodic solitons. The frequency spectrum 
of the radiation shows a remarkable correlation with the 
type of transition the soliton is suffering. 

The indicator (16) with w = 27r/T is positive along 
each curve in Figj2| — moreover, cr(a;) is much greater 
than 7. Therefore, in the case 7 > 0.275, the long-range 
radiation from the soliton is dominated by the first har- 
monic. On the other hand, cr(27r/r) is negative along 
each Tl branch in Figj4] which implies that the first- 
harmonic radiation is short-range when 7 is smaller than 
0.25. (When 7 = 0.10, a decreases from values close to 
—7 to about —67 along the curve; when 7 = 0.20, cr drops 
from zero to nearly —2.57.) In the latter case (7 = 0.10 
and 0.20) the indicator cr(w) becomes positive if one lets 
Lo = 47r/T; hence the radiation is second-harmonic here 
(except for very large T where it is third-harmonic). 

This observation allows us to understand, in qualita- 
tive terms, why the transformation of the soliton with 
7 > 0.275 is different from the scenario realised for very 
weak dampings (7 < 0.25). When 7 is not very small, 
the soliton can dispose the excess energy via two chan- 
nels. Some energy will be lost to the damping in the core 
region, where 1 and the term i^^j) is ~ 7. The 

rest will be sent away in the form of the strong first- 
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harmonic radiation. Due to the availability of this pow- 
erful back-up dissipation channel, the soliton behaves as 
an overdamped system; hence its enhanced stability. On 
the other hand, when 7 is very small, the dominant long- 
range radiation is second harmonic. Since the amplitude 
of the 71-th harmonic radiation scales as e" where e is 
the scale of the amplitude of the oscillation in the core, 
the second-harmonic radiation is much weaker than the 
first-harmonic one, and the radiation channel is effec- 
tively blocked. The soliton shows the volatility of an 
underdamped system; hence its instability and period 
doublings. 

The formation of three-soliton complexes observed in 
the continuation of the oscillatory solitons with 7 — 0.30 
and 7 = 0.35 is also made possible due to the availability 
of strong radiation with the period equal to the period of 
the central soliton. The lateral humps in the complex are 
seeded by the two maxima of the radiation sinusoid which 
are closest to the central soliton. On the other hand, 
the second-harmonic radiation available in the 7 — 0.10 
and 7 = 0.20 cases is unable to seed the lateral solitons 
because it is weak and has the period different from the 
period of the central hump. 

It is interesting to test these qualitative considerations 
using the case of the intermediate damping, 7 = 0.265. 
In this case the first-harmonic indicator cr(27r/T) is posi- 
tive but (7/7 is not much greater than 1 over the most of 
the top Tl branch in Fig[5] (As h grows from 0.35 to 0.71 
to 0.78, the quotient cr /7 on the top branch changes from 
0.5 to to 1.). Accordingly, the first harmonic is strongly 
damped (52 ^ and the long-range radiation is second 
harmonic. As a result, the soliton suffers period-doubling 
bifurcations, as in the case of 7 = 0.10 and 0.20. How- 
ever, as we continue further, the quotient 17/7 grows (to 
2. at /i = 0.82) and the first-harmonic radiation takes 
over from the second-harmonic. Consequently, the solu- 
tion transforms into a three-soliton complex — as in the 
case of 7 = 0.30 and 0.35). 



VII. DISCUSSION AND CONCLUSIONS 
A. Homoclinic explosion 

The numerical continuation of the Tl solution with 
7 < 0.25 shows an unbounded growth of its period. The 
same is true for the T2 solution with 7 < 0.275. These 
observations imply the presence of a homoclinic bifurca- 
tion where the period becomes infinite. We now argue 
that this bifurcation is a source of chaos. 

We use results of Shilnikov [TS] who considered a three- 
dimensional dynamical system with a homoclinic orbit 
connecting a saddle-focus to itself. The saddle-focus is a 
fixed point with a real unstable eigenvalue Ai > and two 
complex-conjugate stable eigenvalues A2 = A3, ReA2,3 < 
0. Assuming that 



Shilnikov proved that at the bifurcation point, the system 
has infinitely many Smale's horseshoes in a neighbour- 
hood of the homoclinic orbit [T5] (see also [IH])- Each 
of these horseshoes contains an invariant Cantor set with 
a countable infinity of (unstable) periodic orbits and an 
uncountable infinity of (unstable) bounded aperiodic or- 
bits — the strange invariant set which manifests itself as 
attractive or transient chaos. Later it was shown that 
finitely many of the horseshoes persist on one side of the 
homoclinic bifurcation |20j . The homoclinic bifurcation 
giving birth to this plethora of periodic and chaotic or- 
bits has been referred to as the "homoclinic explosion" 

The homoclinic explosion may also occur in four- and 
higher-dimensional systems with a homoclinic orbit. The 
necessary conditions for this are formulated in terms of 
the eigenvalues of the linearisation about the fixed point 
connected to itself by the orbit: (a) out of all eigenval- 
ues in the right half of the complex plane, the closest to 
the imaginary axis is a real eigenvalue Ai; (b) in the left 
half of the complex plane, a pair of complex-conjugate 
eigenvalues A2,3 are the closest eigenvalues to the imag- 
inary axis. The homoclinic explosion occurs then if the 
Shilnikov inequality (20) is satisfied pTO]. A carefully 



studied and clearly explained example of homoclinic ex- 
plosion in partial differential equations is in |22| . 

The role of the saddle-focus fixed point in our equation 
([1]) is played by the soliton ?/;_, Eq.([2|. The linear spec- 
trum of this fixed point consists of two discrete eigenval- 
ues, Eqs. ( A7 ) of the Appendix, and a continuum of values 
A = —7 -|- iuj, where |a;| > wq, ujq = Vl — h'^- In the left 
half of the complex plane, the eigenvalue A2 is further 
away from the imaginary axis than the continuous spec- 
trum; hence the dynamics near the saddle-focus will be 
dominated by the single positive eigenvalue Ai and the 
continuum of eigenvalues with negative real part. This is 
not exactly the situation of Shilnikov; however it is simi- 
lar in the sense that the fixed point has a one-dimensional 
unstable manifold and a more-than-one-dimensional sta- 
ble manifold, formed by orbits spiralling into the focus. 
It is natural to expect Shilnikov's result to remain valid 
in this case, provided an analog of the inequality (20) is 
in place — although we do not have any rigorous proof of 
that. The analog of the Shilnikov inequality in our case 
is 



Ai > 7, 



(21) 



with Ai as in (A7) 



The inequality (21) translates into (A^/\/3)Ao > 7, 



where Ao(e) is the positive eigenvalue of the problem 
( A3 ) . One can easily demarcate the region where this 



inequality is valid. From (A2) we have 



Ai > |ReA2,3| 



(20) 



h = 




(22) 
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Equation ( 22 ) together with 



71 



Ao(e) 



(23) 



define the parametric curve h = h^hil) on the (7, /i) 
plane. The curve is parabohc in shape (see the dotted 
curve in the main frame and inset in figj?]); it connects 
the origin to the point 7 = 0, h = l. The inequahty ( 21 ) 



is vahd in the region bounded by this curve and the h- 
axis. This region obviously includes the loci of the homo- 
clinic bifurcations for all 7 < 0.275 (see figjT]). Therefore 
the homoclinic bifurcation is indeed of the homoclinic- 
explosion type in our case. It is this bifurcation that 
is responsible for the appearance of temporally chaotic 
solitons in the parametrically driven damped nonlinear 
Schrodinger equation. 
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B. Attractor charts 

1. Our analysis answers a number of questions raised 
by the attractor chart [12| compiled via direct numer- 
ical simulations. One such question concerns the exis- 
tence of two types of the transformation scenario of the 
parametrically driven soliton, as the driving strength is 
raised for the fixed damping coefficient. Is the trans- 
formation followed by the soliton with large 7 > 0.275 
really different from the transition to chaos suffered by 
its weakly damped counterpart (7 < 0.25), or is this dif- 
ference caused merely by numerical approximations? 

The answer to this question is provided by diagrams 
in Figs. [2] |4j andjsj In the region 7 < 0.25 (represented 
by the values 7 = 0.10 and 0.20, FigQ the increase of 
h results in a sequence of period-doubling bifurcations of 
the periodic soliton. On the contrary, no period-doubling 
bifurcations occur as h is increased for the fixed 7 > 
0.275. (See Fig[2] showing the diagrams for 7 — 0.30 and 
0.35). 

2. Another question is why no periodic solitons are 
observed in numerical simulations with sufficiently large 
h. Is there a well-defined boundary of the domain of exis- 
tence of stable periodic solitons? We have shown that for 
7 > 0.275, the region of existence of periodic solitons on 
the h axis is bounded by the saddle-node bifurcation. No 
such solutions, stable or unstable, exist above this turn- 
ing point. For small 7, 7 < 0.25, the region of existence 
of periodic solitons is also bounded by a turning point; 
however, solitons cease to exist as attractors for lower 
values of h. The boundary of the stability domain here 
is set by a narrow strip of temporally-chaotic solitons. 

These conclusions verify and explain the formation of 
the "desert region" on the (7, h) plane. The positions 
of the saddle-node bifurcation point for 7 > 0.275 and 
the accumulation point of higher-periodic solutions for 
7 < 0.25, are in good agreement with the boundary of 
the desert region detected in simulations [12 . 

3. The third issue clarified by our analysis pertains to 
the shape of the region on the (7, /i)-plane where peri- 



FIG. 7. (Color online) Single-soliton attractor chart. The 
Shilnikov inequality (211 holds inside the region bounded by 
the dotted curve. The inset gives a wider perspective of this 
region. 



odic solitons are observed in numerical simulations. The 
question here is why does the ^1(7) curve bounding this 
region fold on itself for 7 between 0.25 and 0.275. 

This phenomenon is accounted for by the restabilisa- 
tion of the periodic soliton as h is increased for 7 = 0.265 
(Figjsj). We note that the case 7 = 0.265 is intermediate 
between the small- and large-7 transformation scenarios: 
On the one hand, similar to the weakly-damped scenario, 
the soliton undergoes a period-doubling cascade as h is 
raised for this 7 (more precisely, it undergoes two period- 
doubling cascades). On the other, the domain occupied 
by the periodic attractor is bounded by a saddle-node 
bifurcation in this case, as in the strongly-damped situ- 
ation. 

FigjT] summarises our conclusions on single-soliton pe- 
riodic attractors. This diagram is in a good agreement 
with the attractor chart produced using direct numerical 
simulations, Fig[l] 



Appendix A: Eigenvalues of the soliton 

We let -0 = ij)^ {x) + Sip{x, t), where 

Sip = e"**^- [u{x, t) + iv{x, t)] 

is a small perturbation, and linearise Eq.Q in Sip. As- 
suming that 

w(x, t) = Re [e^*-p{x)\ , u(x, t) = Re [e^*g(x)] 

with complex p, q and A, yields an eigenvalue problem 

Lip -f 279 = -\q, {Lo + e)q = \p. (Al) 
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Here we have defined 7 = and A = A/A^, and 

introduced two Sturm-Liouville operators witli familiar 
spectral properties: 

Lo = -d?/dX'^ + 1 - 2sccli^X, 
Li = -d?/dX'^ + 1 -Qscch^X. 



The independent variable X 
e was defined by 



A-x and the parameter 



The two-parameter eigenvalue problem (Al) can be 
reduced [S] to a one-parameter problem 



Lip = -Aq, {Lq + e)q = Ap 



by letting 



A(A + 27) 



(A3) 



(A4) 



and q = {X + 2^)\-'^q. 

The lowest eigenvalue of the operator Lq is zero; hence 
for any e > the operator Lq + e is positive definite and 
the eigenvalue problem ( A3 1 can be cast in the form 



Lip = -A^Mp, 



(A5) 



where M = (Lq + e)"^- The operator in the left-hand 
side is symmetric and the one on the right is symmetric 
and positive definite; hence all eigenvalues — A^ are real 
and the function p{x) can also be considered real. The 
smallest eigenvalue of ( A5 ) can be found as a minimum 
of the Rayleigh quotient: 



Aq = mm 



pe3 {p, Mp) ' 



(A6) 



where (p, g) stands for the scalar product: (p, g) — 
J p{x)q(x)dx. 

The operator Li has a negative eigenvalue, 770 = —3, 
with the eigenfunction yQ = sech^X. Letting p = j/q in 



Eq.(A6), we conclude that the minimum is negative and 
hence the problem (A5| has a negative eigenvalue. We 



denote the associated eigenfunction pq (x) : 
LiPo = -ApMpo- 

Since Li does not have any other negative eigenval- 
ues (the only other discrete eigenvalue of Li is rji = 
0), the problem (A5) does not have any other nega- 
tive eigenvalues either. Indeed, assume there exist — A^ 
(A2) and pi{x) such that Lipi = —A^Mpi. The quadratic 



form (y^Liy) is then negative definite on the subspace 
spanned by po and pi; in particular, it is negative for 
y{x) — CoPq{x) + CiPi{x) where the coefficients are cho- 
sen such that (y, j/o) = 0. This, however, contradicts the 
fact that the form (y, Liy) cannot take negative values 
on the subspace orthogonal to yo. 

Thus the problem ( A3 ) has only one pair of nonzero 
eigenvalues, Aq and — Aq (where Aq is taken to be pos- 
itive). The function Ao(e) has asymptotic behaviours 
Aq ^ (4e)i/2 as e ^ and Aq ^ (3e)i/2 as e — >■ 00: 
for intermediate values of e we have tabulated it numer- 
ically. Using eq.(A4) we recover two eigenvalues of the 



original, two-parameter, problem (Al) 



Al, 



-j±Jj^ + AtAl 



(A7) 



Here Ai > and A2 < -7 < 0. 
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